A SIMPLE AND EFFICIENT NUMERICAL METHOD FOR 
COMPUTING THE DYNAMICS OF ROTATING BOSE-EINSTEIN 
CONDENSATES VIA A ROTATING LAGRANGIAN COORDINATE* 
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Abstract. Wo propose a simple, efficient and accurate numerical method for simulating the 
dynamics of rotating Bose-Einstein condensates (BECs) in a rotational frame with/without a long- 
range dipole-dipolc interaction. We begin with the three-dimensional (3D) Gross-Pitacvskii equation 
(GPE) with an angular momentum rotation term and/or long-range dipole-dipole interaction, state 
the two-dimensional (2D) GPE obtained from the 3D GPE via dimension reduction under anisotropic 
external potential and review some dynamical laws related to the 2D and 3D GPE. By introducing 
a rotating Lagrangian coordinate system, the original GPEs arc rc-formulatcd to GPEs without the 
angular momentum rotation which is replaced by a time-dependent potential in the new coordinate 
system. We then cast the conserved quantities and dynamical laws in the new rotating Lagrangian 
coordinates. Based on the new formulation of the GPE for rotating BECs in the rotating Lagrangian 
coordinates, a time-splitting spectral method is presented for computing the dynamics of rotating 
BECs. The new numerical method is explicit, simple to implement, unconditionally stable and 
very efficient in computation. It is spectral order accurate in space and second-order accurate in 
time, and conserves the mass in the discrete level. Extensive numerical results are reported to 
demonstrate the efficiency and accuracy of the new numerical method. Finally, the numerical method 
is applied to test the dynamical laws of rotating BECs such as the dynamics of condensate width, 
angular momentum expectation and center-of-mass, and to investigate numerically the dynamics and 
interaction of quantized vortex lattices in rotating BECs without/with the long-range dipole-dipolc 
interaction. 
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1. Introduction. Bose-Einstein condensation (BEC), first observed in 1995 [4J 
I17 [ l22 j . has provided a platform to study the macroscopic quantum world. Later, 
with the observation of quantized vortices 018,33,34,36,38,4s], rotating BECs have 
been extensively studied in the laboratory. The occurrence of quantized vortices 
is a hallmark of the superfluid nature of Bose-Einstein condensates. In addition, 
condensation of bosonic atoms and molecules with significant dipolc moments whose 
interaction is both nonlocal and anisotropic has recently been achieved experimentally 
in trapped 52 Cr and 164 Dy gases [HHHIiniEIlEllESlIlS] ■ 

At temperatures T much smaller than the critical temperature T c , the proper- 
ties of a BEC in a rotating frame with long-range dipole-dipole interaction are well 
described by the macroscopic complex-valued wave function </' = t/'(x. t), whose evo- 
lution is governed by the three-dimensional (3D) Gross-Pitaevskii equation (GPE) 
in dimcnsionless units with angular momentum rotation term and long-range dipolc- 
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dipole interaction OQIHIiniHniliniHZlISS] : 



(1.1) 



id t tp{x.,t) 



--V 2 + V(x) + k#| 2 + A ([/ dip * H 2 ) - QL 2 



where t denotes time, x = (x, y, z) T G R 3 is the Cartesian coordinate vector, V(x) 
is a given real-valued external trapping potential which is determined by the type of 
system under investigation and K = is a dimcnsionless constant describing the 

strength of the short-range two-body interaction (positive for repulsive interaction, 
and resp. negative for attractive interaction) of a condensate consisting of N particles 
with s-wave scattering length a s and a dimcnsionless length unit x s [9l[T9]. Further- 
more £1 G R is a given dimcnsionless constant representing the angular velocity, and 

A = 



3K 2 x s 



is a dimensionless constant describing the strength of the long-range 
dipolc-dipolc interaction with m the mass of a particle, fio the vacuum permeability, 
/^dip the permanent magnetic dipole moment (e.g. [idip = 6/^ B for 52 C r with /x B being 
the Bohr magneton) and ft the Planck constant [HlUniHD] ■ I n addition, L z is the 
dimensionless z-component of the angular momentum rotation defined as 



(1.2) 



-i(xd y - yd x ), 



and C/dip is the dimcnsionless long-range dipolc-dipolc interaction potential defined as 



(1.3) £/ di p(x) 



4-7T X 1 



3(x • n) 



4-7T X ' 



[l - 3cos 2 (i9)] 



x G 



with n = (711, ri2, "3) T <G R 3 a given unit vector, i.e. |n| = y/ n\ + n 2 . + n\ = 1, 
representing the dipole axis (or dipole moment) and 7? = $ n (x) the angle between the 
dipole axis n and the vector x [9}; 19| . The wave function is normalized to 



(1.4) 



IH| 2 := 



/ h/>(x,*)| 2 dx = l. 



In addition, similar to [9j|T9] , the above GPE (jl.ip can be re- formulated as the fol- 
lowing Gross-Pitaevskii-Poisson system [6ll9lll9] 



(1.5) id t ^t) 



-V 2 + V(x) + (k- A) I VI 2 - 3A</?(x,i) - &L 2 



V>(x,i), 



(1.6) <p(x,t) = d an u(x,t), -V 2 w(x,i) = ^(x,*)! 2 with lim u(x,i)=0, 



where <9 n = n • V and 9, 
1 



(1.7) u(x,t) 



47T X 



|x|— ^00 

9 n (9 n ). From (|1.6|) . it is easy to see that 
1 



4-7T X - X' 



-|V(x',*)|W, 



x G 



t > 0. 



In some physical experiments of rotating BEC, the external trap is strongly con- 
fined in the z-dircction, i.e. 



(1.8) 



V(x) = V 2 (x,y) + 



2e 4 ' 



x € 



with < £ <C 1 a given dimensionless parameter [BJ, resulting in a pancake-shaped 
BEC. Similar to the case of a non- rotating BEC, formally the GPE ([TT]) or (TT3)) - (|TTd] 
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can effectively be approximated by a two-dimensional (2D) GPE as [51 [1T71[T^] : 



(1.9) id t Tp(x±,t) 



L, _., , k + A(3n§ - 1) . , l2 3A 



(1.10) p = p(x ±> f) = (dnxnx -niVi)«(x ±J t), x ± = (x,zj) T eM 2 , i>0, 

where V± = (d x ,d y ) T , \7\ = d xx + d yy , n ± = (n 1 ,n 2 ) T ', d n± = n ± • V_l, d n 
dn ± (d n± ) and 

(1.11) u(x ± ,t) = G £ *\ij\ 2 , G e (x ± ) 



-s 2 /2 



ids, xj_ € 



The above problem (|1.9[) - (|1.10[) with (JTTTTJ) is usually called surface adiabatic model 
(SAM) for a rotating BEC with dipole-dipole interaction in 2D. Furthermore, taking 
£ — s- 0+ in (|l.lll) , we obtain 



(1.12) 



1 



G°(xi)> x ± g 



2tt|x_l 

This, together with (jl.lip . implies that when e — > + 
1 



(1.13) u(xj.,t) 



2tt xi 



(-Vi) 1/2 u = |0| 2 with lim u(x_L,i) = 0. 

|xj_|->oo 



The problem (|1.9I) - (|1.10I) with (|1.13[) is usually called surface density model (SDM) 
for a rotating BEC with dipole-dipole interaction in 2D. Note that even for the SDM 
we retain the e-dependence in (|1.9[) . 

In fact, the GPE (JTTJ or (JOJ in 3D and the SAM or SDM in 2D can be written 
in a unified way in d-dimensions (d = 2 or 3) with x = (x, y) T when d = 2 and 
x = (x, y, z) T when d = 3: 



(i.i4) ia^(x,t) 



iv 2 + y(x) + /3|v^| 2 + w (x, t) - fii 2 



(1.15) p(x,i) =L„«(x,*), u(x,i) = G* IV'I 2 , xGM d , i > 0, 
where V(x) = V2(x) when d = 2 and 



(1.16) /? 



K+A(3n§-1) 
k — A, 



V = 



-3A/2, 
-3A, 



in = 



On 
8 n 



n|V 2 , d = 2, 
d = 3, 



1/2tt|x|, 
(1.17) G(x) = <( G £ (x), 
1/4tt|x|, 



G(0 



i 



27/iW ds ' d = 2&SAM, 



d = 2 & SDM, 
d = 2 
d = 3, 



where /(£) denotes the Fourier transform of the function /(x) for x, £ <G M d . For 
studying the dynamics of a rotating BEC, the following initial condition is used: 



(1.18) V(x,0)=Vo(x), 



x g 



with ||V>o|| 



|Vo(x)| 2 dx=l. 
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Wc remark here that in most BEC experiments, the following dimcnsionlcss harmonic 
potential is used 



where j x > 0, "f y > and 7 Z > are dimensionless constants proportional to the 
trapping frequencies in x-, y- and z-direction, respectively. 

Recently, many numerical and theoretical studies have been done on rotating 
(dipolar) BECs. There have been many numerical methods proposed to study the 
dynamics of non-rotating BECs, i.e. when fl = and rj = [3|flT | [14 l [20 l [28 l [37 1 l43] . 
Among them, the time-splitting sine/Fourier pseudospectral method is one of the most 
successful methods. Compared to other methods, the time-splitting pseudospectral 
method has spectral accuracy in space and is easy to implement. In addition, as was 
shown in [5], this method can also be easily generalized to simulate the dynamics of 
dipolar BECs when rj ^ 0. However, in rotating condensates, i.e. when ft ^ 0, we can 
not directly apply the time-splitting pseudospectral method proposed in [14) to study 
their dynamics due to the appearance of angular rotational term. So far, there have 
been several methods introduced to solve the CPE with an angular momentum term. 
For example, a pseudospectral type method was proposed in |10j by reformulating the 
problem in the two-dimensional polar coordinates (r, 0) or three-dimensional cylindri- 
cal coordinates (r, 9, z). The method is of second-order or fourth-order in the radial 
direction and spectral accuracy in other directions. A time-splitting alternating di- 
rection implicit method was proposed in [13] , where the authors decouple the angular 
terms into two parts and apply the Fourier transform in each direction. Furthermore, 
a generalized Laguerre- Fourier- Hermite pseudospectral method was presented in [12] . 
These methods have higher spatial accuracy compared to those in J3[[8j[28] and are 
also valid in dissipative variants of the GPE (|1.1[) . cf. [33]. On the other hand, the 
implementation of these methods can become quite involved. The aim of this paper 
is to propose a simple and efficient numerical method to solve the GPE with angular 
momentum rotation term which may include a dipolar interaction term. One novel 
idea in this method consists in the use of rotating Lagrangian coordinates as in [5] in 
which the angular momentum rotation term vanishes. Hence, we can easily apply the 
methods for non-rotating BECs in [JJ] to solve the rotating case. 

This paper is organized as follows. In Section [21 we present the dynamical laws 
of rotating dipolar BECs based on the GPE (|iT4l) - ([i~li?|) . Then in Section H we 
introduce a coordinate transformation and cast the GPE and its dynamical quantities 
in the new coordinate system. Numerical methods are proposed in Section [3] to 
discretize the GPE for both the two-dimensional and three-dimensional cases. In 
Section we report on the accuracy of our methods and present some numerical 
results. We make some concluding remarks in Section [B] 

2. Dynamical properties. In this section, we analytically study the dynamics 
of rotating dipolar BECs. We present dynamical laws, including the conservation of 
angular momentum expectation, the dynamics of condensate widths and the dynamics 
of the center of mass. In the following we omit the proofs for brevity; they are similar 
to the ones in [10l[14] . 

2.1. Conservation of mass and energy. The GPE in (|1.14[) - (|1.18[) has two 

important invariants: the mass (or normalization) of the wave function, which is 



(1.19) 




d = 2, 
d = 3, 
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defined as 

(2.1) N(t) := !|^(-,t)|| 2 := / |^(x, t)| 2 dx = f |^(x, 0)| 2 dx = 1, t > 0, 

JR d JR d 

and the energy per particle 

i|vy,| 2 + y(x)^ 2 + |m 4 + |^| 2 - nrL z ip 



dx 



E(t) :=E(i>(;t)) 

(2.2) = £(</;(• ,0))=£(V>o), i>0, 

where /* denotes the conjugate of the complex-valued function /. Stationary states, 
corresponding to critical points of the energy defined in (|2.2[) , play an important role 
in the study of rotating dipolar BECs. Usually, to find stationary states s (x), one 
can use the ansatz 

(2.3) 4>(x,t) =e- Vst s (x), xeR d , t>0, 

where /j, s <G M is the chemical potential. Substituting (|2.3p into (|1.14[) yields the 
nonlinear eigenvalue problem 

(2.4) Ms0 s (x) = 



-iv 2 + V{x) + /?|0 S | 2 + Wa - QL Z 



(2.5) ip s (x) = L n u s {x), m s (x) = G* \cj) s \ 2 7 xe 
under the constraint 



c^ s (x), xe 

3)'' 



(2.6) ||</> s || 2 = / |^(x)|^x = l. 

JR d 

Thus, by solving the constrained nonlinear eigenvalue problem (|2.4p ~ (|2.6p . one can 
find the stationary states of rotating dipolar BECs. In physics literature, the station- 
ary state with the lowest energy is called ground state, while those with larger energy 
are called excited states. 

2.2. Conservation of angular momentum expectation. The angular mo- 
mentum expectation of a condensate is defined as |10U14) 

(2.7) (L z )(t)= [ iP*{x,t)L z iP(x,t)dx, t>0. 

JR d 

This quantity is often used to measure the vortex flux. The following lemma describes 
the dynamics of angular momentum expectation in rotating dipolar BECs. 

Lemma 2.1. Suppose that ip(x,t) solves the GPE \1.1^ - U.18\) with V(x) chosen 
as the harmonic potential \1.19\) . Then we have 

(2-8) &M = ( 7 2 - 7 2 )/ xy\4>\ 2 dx-ri[ \^\ 2 [{xd v -yd x ) V ]dx, t > 0. 

al JR d JR d 

Furthermore, if the following two conditions are satisfied: (i) j x = j y and (ii) rj = 
for any given initial data ipo in HUM) or n — (0,0, 1) T when i[>o satisfies ipo(x) = 
f(r)e lmB in 2D and V'o(x) = f( r , z)e lm9 in 3D with m <E 1 and f a function of r in 
2D or (r,z) in 3D, then the angular momentum expectation is conserved, i.e. 

(2.9) (L z )(t) = (L z )(0), t>0. 

That is, in a radially symmetric trap in 2D or a cylindrically symmetric trap in 
3D, the angular momentum expectation is conserved when either there is no dipo- 
lar interaction for any initial data or the dipole axis is parallel to the z-axis with a 
radially /cylindrically symmetric or central vortex-type initial data. 
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2.3. Dynamics of condensate width. The condensate width of a BEC in a- 
dircction (where a = x,y,z or r~ \/ x 2 + y 2 ) is defined by 

(2.10) a a {t) = y/djfj, t>0, where 8 a (t) = a 2 \ip(x, i)| 2 dx. 



In particular, when d = 2, we have the following lemma for its dynamics |10j : 

Lemma 2.2. Consider two-dimensional BECs with radially symmetric harmonic 
trap U.19\) , i.e. d = 2 and 7:r = j y := 7r . If r\ = 0, then for any initial datum ipofa) 
in | fPgp . it holds for t > 

(2.11) S r (t) = + gi^PJ [1 _ C0S (2 7r t)] + 8™ cos(2 7r t) + S in(2 7r t), 

7r 2 7r 

w/iere <5 r (i) := 6 x (t) + <^ 0) := 4(0) + 8 y (0) and 8^ := 4(0) + 4,(0). Further- 

more, if the initial condition "00 ( x ) *s radially symmetric, we have 

1 

— < 
2 

r(l) 



= a tf (t) - -8 r (t) 



(2.12) = [i _ cos (2 7lt) ] + ^) C0S (2 7x t) + ^ sin(2 Tl t), t > 0. 

Thus, in this case the condensate widths o~ x (t) and o~ y (t) are periodic functions with 
frequency doubling trapping frequency. 

2.4. Dynamics of center of mass. We define the center of mass of a conden- 
sate at any time t by 

(2.13) x c (t) = / x|?/>(x,t)| 2 dx, <>0. 

The following lemma describes the dynamics of the center of mass. 

Lemma 2.3. Suppose that ip(x,t) solves the GPE \1.14\ j- M.18}) with V(x) chosen 
as the harmonic potential (1.19\) . Then for any given initial data ipo, the dynamics of 
the center of mass are governed by the following second-order ODEs: 

(2.14) x c (t) - 2f2Jx c (i) + (A + Q 2 J 2 )x c (t) = 0, t>0, 

(2.15) x c (0)=x£°) := f x|^ | 2 rfx, 

(2.16) x c (0) = xW := f Im(i/;*V^o)rfx-f7Jx( 0) , 

where Im(/) denotes the imaginary part of the function f and the matrices 

for d = 2, 



or 







-1 
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for d = 3. 



Lemma l2.3l shows that the dynamics of the center of mass depends on the trapping 
frequencies and the angular velocity, but it is independent of the interaction strength 
constants /3 and rj in (|1.14[) . For analytical solutions to the second-order ODEs (|2.14[) - 
(|2~Td| . we refer to [ST]. 
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2.5. An analytical solution under special initial data. From Lemma [ 
we can construct an analytical solution to the GPE (|1.14p — 8[) when the initial data 
is chosen as a stationary state with its center shifted. 

Lemma 2.4. Suppose V(x) in |7T7ZP * s chosen as the harmonic potential h!.19\ ) 
and the initial condition ?/>o( x ) Mi hl.l&j) is chosen as 

(2.17) Vo(x) =0 s (x-x°), xeM d , 

where x° € K d is a given point and </> s (x) is a stationary state defined in {2.4\ )- [2~6\ ) 
with chemical potential fi s , then the exact solution of \1. lJ$ - ll.l8)) can be constructed 
as 

(2.18) = s (x-x c (t))e-^*e ra(x ^ ) , x e M d , i > 0, 
where x c (i) satisfies the ODE {2.1$ with 

(2.19) x c (0) = x°, x c (0) = -ftJx , 
and w(x,t) is linear in x, i.e. 

iy(x,i) =c(t) -x + g(t), c(t) = (ci(i), . . . , c d (t)) T , x e R d , t>0 

for some functions c(t) and g(t). Thus, up to phase shifts, ip remains a stationary 
state with shifted center at all times. 

3. GPE under a rotating Lagrangian coordinate. In this section, we first 
introduce a coordinate transformation and derive the GPE in transformed coordi- 
nates. Then we reformulate the dynamical quantities studied in Section [5] in the new 
coordinate system. 

3.1. A rotating Lagrangian coordinate transformation. For any time t > 
0, let A(t) be an orthogonal rotational matrix defined as 

(3.1) A(t) = ( C0S % Sh « V if d = 2, 
v ' V — 8in(S2t) cos(iM) J 

and 

/ cos(fii) Bin(fii) \ 

(3.2) A(t) = - sin(fii) cos(ftt) , if d = 3. 

V o 01/ 

It is easy to verify that j4 _1 (t) = A T {i) for any t > and A(0) = / with / the identity 
matrix. For any t > 0, we introduce the rotating Lagrangian coordinates x as [5j[23,27 

(3.3) x = A _1 (t)x = ,4 T (t)x & x = A(i)x, xeM d , 
and denote the wave function in the new coordinates as <fi := </>(x, t) 

(3.4) <f>(5c,t) :=ip(-x.,t) =ip(A(t)x,t) , x <G M. d , t>0. 

In fact, here we refer the Cartesian coordinates x as the Eulerian coordinates and 
Fig. 13. II depicts the geometrical relation between the Eulerian coordinates x and the 
rotating Lagrangian coordinates x for any fixed t > 0. 
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V ^ 


y 

X 

, ,T 

, \ 'sit 

' — 1 > 




x a: 



Fig. 3.1: Cartesian (or Eulerian) coordinates (x,y) (solid) and rotating Lagrangian 
coordinates (x, y) (dashed) in 2D for any fixed t > 0. 



Using the chain rule, we obtain the derivatives 



a^(x,i) = dtH^t) + W(x,t) • (4(t)x) = atV(x,t) - n(atf B - yd x )^(x,t), 
V0(x,<) = J 4~ 1 (<)VV'(x,i), V 2 0(x,t) = V 2 </>(x,i). 

Substituting them into (|1.14|l - (|1.18|) gives the following d-dimensional GPE in the 
rotating Lagrangian coordinates: 



(3.5) idt(f>(x,t) = 

(3.6) (f(x,t) = L m{t) u(x,t), u(x,t) = G * |0| 2 , 
where G is defined in (|1.17|l and 

(3.7) W(x,t) = V{A(t)x), xeK d , 



x G 



x,t), 



t > 0, 



t > 0, 



(3.8) 



*(t) 



mj_(t)mj_(t) 
m(i)m(i) 5 



ti 2 V 2 , d=2, 



d = 3, 



f > 0, 



with m(t) G R 3 and mj_ (t) G M 2 defined as 
mi(t) 

(3.9) m(i) = I m 2 (*) I := (t)n = I m sin(Qi) + n 2 cos(fii) 



ni cos(f2t) — ri2 sin(Ot) 



t > 0, 



m 3 (t) 



and mj_(f) = (mi (t), m2(t)) T , respectively. The initial data transforms as 



(3.10) 



(x,0) = V'CxjO) = Vo(x) := <fo(x) = </>o(x), 



x = x G 



We remark here again that if V(x) in p,14[) is a harmonic potential as defined in 
(fj~T9l) , then the potential W(x,t) in ([3~5]) has the form 



W / (x,t) = i 



wi(z 2 + y 2 ) + w 2 
wi(x 2 +y 2 ) + ui 2 



{x 2 - y 2 ) cos(2fti) + 2xy sin(2Qt) 
(x 2 - y 2 ) cos(2fti) + 2xy sin(2fit) 



d = 2, 
+ 27^, d = 3, 
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where uji = 7 2 + 7 2 and u> 2 = 7 2 ~ 7y- It is easy to see that when j x — i y := 7 r , 
i.e. radially and cylindrically symmetric harmonic trap in 2D and 3D, respectively, 
we have u>i = 2j 2 and 102 = and thus the potential W(x, £) = V^(x) becomes 
time- independent . 

In contrast to (| 1 . 14|) . the GPE (j33|) does not have an angular momentum ro- 
tation term, which enables us to develop simple and efficient numerical methods for 
simulating the dynamics of rotating dipolar BEC in Section |4] 

3.2. Dynamical quantities. In the above, we introduced rotating Lagrangian 
coordinates and cast the GPE in the new coordinate system. Next we consider the 
dynamical laws in terms of the new wave function </>(x, t) . 

Mass and energy. In rotating Lagrangian coordinates, the conservation of mass 
(ED) yields 



(3.11) ||V(-,*)II 2 := / |^(x,i)| 2 dx= / \4>(x,t)\ 2 ( &=U(;t)\\ 2 = l, t>0 



dx 



The energy defined in (|2.2I) becomes 



i|V^| 2 + M/(x,t)|0| 2 + ^| 4 + ^|0| 2 



E(c/>(;t)) -- 

(3.12) -f f [9 r ^(x,r) + ^(a T i m(T) ) W (x,r)l|0| 2 drdx = ^(0(-,O)), t > 0, 
Jm d Jo 1 1 J 



where u is given in (|3.6[) . Specifically, it holds 



Angular momentum expectation. The angular momentum expectation in the 
new coordinates becomes 

(L z )(t) = -i f V*(x,t)(^-y5 x )V(x,t)dx 

(3.13) = — i / 0*(x, t)(xd y — y9j)(/>(x, t) c£x, t > 0, 

which has the same form as ()2.7|) in the new coordinates of x G K d and the wave 
function </>(x, t). Indeed, if we denote L z as the ^-component of the angular momentum 
in the rotating Lagrangian coordinates, we have L z = —i(xdy — ydx) = —i(xd y — 
yd x ) = L z , i.e. the coordinate transform does not change the angular momentum in 
z-direction. In addition, noticing that for any t > it holds 0(x,t) = ip(x,t) and 
\A(t)\ = 1 for any t > immediately yields (|3.13j) . 

Condensate width. After the coordinate transform, it holds 

(3.14) S r (t)= [ (x 2 +y 2 M 2 dx= f (x 2 +y 2 M 2 dx = 5 x (t)+S v (t), 

(3.15) S z (t)= I z 2 \il>\ 2 dx= [ ^|0| 2 dX = 5 s (t), 
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for any t > 0. 

Center of mass. The center of mass in rotating Lagrangian coordinates is defined 
as 

(3.16) x c (t) = / x|0(x,i)| 2 dx, t>0. 

Since det(A(t)) = 1 for any t > 0, it holds that x c (i) = A(t)x c (t) for any time t > 0. 
In rotating Lagrangian coordinates, we have the following analogue of Lemma l2.4t 

Lemma 3.1. Suppose V(x) in J 1. 11$ is chosen as the harmonic potential M.19\ ) 
and the initial condition </>o(x) in i3.10\) is chosen as 



(3.17) 0o(X) = ^(x-x o ), xe 



where x° is a given point in R d and 4> s (x) is a stationary state defined in l[2.4\ )- ^276\ 
with chemical potential fx s , then the exact solution of Q3.5P — (|3.6[> is of the form 

(3.18) <j>(x,t) = s (x-x c (i))e- Vs *e™ ( ^ t) , f>0, 
where x c (t) satisfies the second-order ODEs: 

(3.19) x c (t) + A T (i)AA(t) x c (i) = 0, i > 0, 

(3.20) x c (0) = x°, x c (0) = 0, 



wif/i £/ie matrix A defined in Lemma \2.3\ and w(x,t) is linear in x, i.e. 

w(x,t) = c(t) -x + g(t), c(t) = (ci(i), . . . ,c d {t)) T , x e R d , i > 0. 

We have seen that the form of the transformation matrix A(t) in (|3.2j) is such 
that the coordinate transformation does not affect the quantities in z-direction, e.g. 
(L z ), <r z (t) and z c (t). 

4. Numerical methods. To study the dynamics of rotating dipolar BECs, in 
this section we propose a simple and efficient numerical method for discretizing the 
GPE p.5[) - p.l0[) in rotating Lagrangian coordinates. The detailed discretizations for 
both the 2D and 3D GPEs are presented. Here we assume £7 ^ 0, and for fi = 0, we 
refer to [7||1QI1II6] . 

In practical computations, we first truncate the whole space problem (|3.5p - (|3.10l) 
to a bounded computational domain T> C R d and consider 

(4.1) id t (f>(x,t) = -^V 2 (j) + W{x,t)<P + f3\<f>\ 2 <t> + ricp<f>, x e V, t > 0, 

(4.2) tp(x,t) = L m{t) u(x,t), u(x,t)= G(x-y)p(y,t)dy, x e V, t > 0; 



0.1 -\ F' !)|2> L e Jl„ feu*. 



where 

0, otherwise 
The initial condition is given by 

(4.3) 0(x,O) = <Mx), xeP. 
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The boundary condition to (|4.1[) will be chosen based on the kernel function G de- 
fined in (|1.17p . Due to the convolution in (|4.2[) . it is natural to consider using the 
Fourier transform to compute it(x, t). However, from (|1.17[) and (|3.1ip . we know 

that lim^o = 00 an d |</>| 2 (£ = 0) ^ 0. As noted for simulating dipolar BECs 
in 3D [9l|15II39|. there is a numerical locking phenomena, i.e. numerical errors will 
be bounded below no matter how small the mesh size is, when one uses the Fourier 
transform to evaluate u(5c,t) and/or ip(x, t) numerically in (|4.2[) . As noticed in 011], 
the second (integral) equation in (|4.2j) can be reformulated into the Poisson equation 
(II. Gp and square-root-Poisson equation (|1.13|) for 3D and 2D SDM model, respec- 
tively. With these PDE formulations for u(x, t), we can truncate them on the domain 
T> and solve them numerically via spectral method with sine basis functions instead of 
Fourier basis functions and thus we can avoid using the 0- modes [3] . Thus in 3D and 
2D SDM model, we choose the homogeneous Dirichlet boundary conditions to (|4.1[) . 
Of course, for the 2D SAM model, one has to use the Fourier transform to compute 
u(x,t), thus we take the periodic boundary conditions to (|4.ip . 

The computational domain V C M d is chosen as T> = [a, b] x [c, d] if d = 2 and 
T> = [a, b] x [c, d] x [e, /] if d = 3. Due to the confinement of the external potential, 
the wave function decays exponentially fast as |x| — > 00. Thus if we choose T> to be 
sufficiently large, the error from the domain truncation can be neglected. As long as 
we solve 0(x, t) in the bounded computational domain T>, we obtain a corresponding 
solution ip(yi,t) in the domain A(t)T>. As shown in Fig. 14.11 for the example of a 2D 
domain, although the domains A(t)T> for t > 0, are in general different for different 
time t, they share a common disk which is bounded by the inner green solid circle in 
Fig. 14.11 Thus, the value of tp{'x.,i) inside the vertical maximal square (the magenta 
area) which lies fully within the inner disk can be calculated easily by interpolation. 




A(t)V 



Fig. 4.1: The bounded computational domain T> in rotating Lagrangian coordinates x 
(left) and the corresponding domain A(t)T> in Cartesian (or Eulerian) coordinates x 
(right) when f2 = 0.5 at different times: t = (black solid), t = j (cyan dashed), t = ? 
(red dotted) and t = ^ (blue dash-dotted). The two green solid circles determine 
two disks which are the union (inner circle) and the intersection of all domains A(t)T> 
for t > 0, respectively. The magenta area is the vertical maximal square inside the 
inner circle. 



4.1. Time-splitting method. Next, let us introduce a time-splitting method 
to discretize (|4.ip - (|4.3p . We choose a time-step size At > and define the time 
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sequence as t n = nAt for n£N. Then from £ = t n to £ = £„+i, we numerically solve 
the GPE (|4.ip in two steps. First we solve 

(4.4) id t cj>(Z,t) = --V 2 <l>(x,t), xel?, t n <t<t n+1 
for a time step of length At, and then we solve 

(4.5) idt<f>{x,t) = [W(x,t)+f3\^\ 2 + w ] 0(x,£), x€P, t n <t<t n+1> 

(4.6) <p(x.,t) = L m(t) u(x,£), u(x,t)= G(x-y)p{y,t)dy, 

Jm d 

for the same time step. 

Equation (|4.4[) can be discretized in space by sine or Fourier pseudospectral 
methods and then integrated exactly in time. If homogeneous Dirichlet boundary 
conditions are used, then we choose the sine pseudospectral method to discretize it; 
otherwise, the Fourier pseudospectral method is used if the boundary conditions are 
periodic. For more details, see e.g. [9|fl4]. 

On the other hand, we notice that on each time interval [£„,£„ + i], the problem 
(14. 5p (|4.6j) leaves |0(x,£)| and hence u(x,£) invariant, i.e. |0(x, £)| = |0(x,£ n )| and 
w(x, £) = u(x, £ n ) for all times t n < t < £ n +i. Thus, for t G [£„,£„ + i], Eq. (|4.5[) 
reduces to 

(4.7) id t <f>{x,t) = [W(x,t) + P\<f>{x,t n )\ 2 + 77 (L m(t)M (x,£„))] <f>{x,t), xeP. 
Integrating f|4. 7[) in time gives the solution 



|2/ 



(4.8) 0(x,£) = 0(x,£„)cxp -i (£|0(X,£ n )| 2 (£-£ n )+r7$(x,£) + y W(x,t)<Zt 
for x £ D and £ € [£„,£„ + i], where the function $(x,£) is defined by 

(4.9) $(x,£) = ^ [L m{r] u(x,t n )]dT = ^ L m{T} dr^j u(x,t n ). 
Plugging (|3l?l) and (f3T8|) into (|3l9]) . we get 

(4.10) $(x,£) =L d (£)u(x,£„), £n<£<£n+l, 

where 



L d (t) 
with 



[ll\t) - lf{t)]d iS: + [If it) - lf(t)]dyy + ; e 12 (£)^, d = 2, 

Wdi* + lfit)dyy + lf{t)d~ Z - Z + lf{t)d. iy + Ifim-Z + lf{t)d m d = 3, 



jT m\{r)dT = [nf cos 2 (fir) + n 2 sin 2 (Sir) - nin 2 sin(20r)] c£r 

^I±^(£ - £ n ) + ^i^t [sin(2ft£) - sin(2tt£ n )] + ^ [cos(2r!£) - cos(2tt£„)] 



/ 22 (t) = jf m|(r)dr = [n 2 , cos 2 (fir) + n 2 sin 2 (Sir) + nin 2 sin(20r)] dr 



7 ^l±^(t - £„) - [sin(2n£) - sin(2ft£ n )] - ^ [cos(2r!£) - cos(2tt£„)] 
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2 mi(r)m2(r)dr = [(n? - n|) sin(2Qr) + 2nin 2 cos(2Qt)] 

2 2 

[cos(2fti„) - cos(2ftt)] + -^—=- [sin(2Qt) - sin(20i„)] , 



;13 



ig (i) = 2n 3 j mi(r)dT = 2n 3 y [ni cos(f2r) — n 2 sin(fir)] dr 

2^3 

= [m [sin(Oi) — sin(fit„)] + n 2 [cos(f2t) — cos(Oi n )]] , 

Z 23 (i) = 2n 3 m 2 (r)dr = 2n 3 [m SHi(Ot) + n 2 cos(Qt)] dr 
2n 3 

[ri2 [sin(fif) — sin(fit„)] — ni [cos(fit) — cos(f2i„)]] , 



i 

2 



l e(t)= nidr = ni{t-t n ), t n <t<t 



n+l- 



In Section 321 we will discuss in detail the approximations to $(x,i) in (|4.10p . In 
addition, we remark here again that if V(x) in (|1.14p is a harmonic potential as defined 
in (|1.19p . then the definite integral in (|4.8p can be calculated analytically as 

where 

1 /"* 

H (x, t) = - y cj 2 [(x 2 - j/ 2 ) cos(2^t) + 2xy sin(20r)] dr 

= ^ [(x 2 - y 2 ) [sin(2fit) - sin(2fii„)] - 2xy [cos(20i) - cos(2Qt n )]l . 

Of course, for general external potential V(x) in (| 1 . 14|) . the integral of W(Sc, r) in (14.81) 
might not be found analytically. In this situation, we can simply adopt a numerical 
quadrature to approximate it, e.g. the Simpson's rule can be used as 



J W(5t,T)dT 



W{Z, tn) + iW(x, + W{St, t) 



We remark here that, in practice, we always use the second-order Strang splitting 
method [42] to combine the two steps in f|4.4[) and (|4.5p - (|4.6p . That is, from time 
t = t n to t = t n+ i, we (i) evolve (|4.4p for half time step At/2 with initial data given 
at t = t n \ (ii) evolve (|4.5p ~- (|4.6j) ) for one step At starting with the new data; and (hi) 
evolve (|4.4p for half time step At/2 again with the newer data. For a more general 
discussion of the splitting method, we refer the reader to [7JIH1I2I] ■ 

4.2. Computation of $(x, t). In this section, we present approximations to the 
function $(x, t) in (|4.10p . From the discussion in the previous subsection, we need 
only show how to discretize u(x, t n ) in (|4.2p and its second-order derivatives in (|4.10p . 

4.2.1. Surface adiabatic model in 2D. In this case, the function u(x, t n ) in 
(|4.9p is given by 



(4.11) u(x,t n )= [ G(x-y)p(y,i n )dy, xeD. 

JR 2 
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with the kernel function G defined in the second line of (|1.17l) . To approximate it, we 
consider a 2D box T> with periodic boundary conditions. 

Let M and K be two even positive integers. Then we make the (approximate) 
ansatz 

M/2-1 K/2-1 

(4.12) u(%t n )= ]T Yl ^(^)e^ (£ - a) e^ (27 - c) , x=(£,y)eX>, 

p=-M/2 q=-K/2 

where u^ q {t n ) is the Fourier coefficient of u(5c,t n ) corresponding to the frequencies 
(p, q) and 

The index set Smk is defined as 

f , m M M K K 

S MK = \(p,q) \ - T <P<y-l, -y<9<y-l 

We approximate the convolution in (|4. by a discrete convolution and take its 
discrete Fourier transform to obtain 

(4.13) uf pq (t n ) = G(^> 2 ) ■ (\H 2 ) f Pg , (p, l) G <W, 

where (|^ ra | 2 )p 9 is the Fourier coefficient corresponding to the frequencies (p, q) of the 
function |0(x, t n )| 2 , and G{v^,v q ) are given by (see details in (|1.17[l ) 

(4.14) G(„ p1 „,) = — 2 J ^ K)2 + K)2 + g2 ^ (p, 9) e 

Since the integrand in (|4.14p decays exponentially fast, in practice we can first truncate 
it to an interval [si,s 2 ] with |si|,S2 > sufficiently large and then evaluate the 
truncated integral by using quadrature rules, e.g. composite Simpson's or trapezoidal 
quadrature rule. 

Combining (|4.10j) . (|4.12|) and (|4.13j) . we obtain an approximation of $(x, t) in the 
solution (|4.8j) via the Fourier spectral method as 



M/2-1 K/2-1 



(4.15) *(x,t)= E [L{viy q ,t)G{ul,ul)-{\^\ 



\2)f 
/pq 

p=-M/2 q=-K/2 '" 



ii/p(x-a) J,v\(y-c) 



for time t n < t < t n +i, where the function L(£i,£2,i) is defined as 

= ~ [WHt) ~ «"(*)) £ + (lf(t) - !»(*)) ?! + Z?(t)&6] • 

4.2.2. Surface density model in 2D. In this case, the function u(x,t n ) in 
(|4.9|) also satisfies the square-root-Poisson equation in (|1.13[) which can be truncated 
on the computational domain T> with homogeneous Dirichlet boundary conditions as 

(4.16) (-V 2 ) 1/2 «(x>*n) = |0(x,t„)| 2 , xeP; u(x,t„)|au = 0. 
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The above problem can be discrctized by using a sine pseudospectral method in which 
the 0-modes are avoided. Letting M, K £ N, we denote the index set 

Tmk = {(p,q) 1 1 <P < M- 1, l<q<K-l}, 

and define the functions 

U Pt g(x) = sin(/itp(x - a))sin(// 2 (j/- c)), (p, q) £ Tn/if , x = (x, j/) £ £>, 

where 

(4.17) ^1 = -^-, ^2 = (p,q)&T M i<- 

F o — a ' a — c 

Assume that 

(4.18) u(x,t n ) = E Z "w^^W*)' * e P ' 

where u pq (t n ) is the sine transform of u(x, i„) at frequencies (p, q). Substituting (|4.18[) 
into (|4.16[) and taking sine transform on both sides, we obtain 

(4.19) u pq {t n ) - Uy ' ^ , (p, q) £ Tmk, 

yW + W 

where (l^™! 2 )^ is the sine transform of |</)(x, i„)| 2 at frequencies (p,q). 

Combining (|4.18[) . (|4. 19[) and f|4. 10[) . we obtain an approximation of $(x, t) in the 
solution (|4.8[) via sine spectral method as 

(4.20) $(x,t) = E E / = [i(^,^,0^, 9 (x) + ?l 2 W^(x)] , 
where the functions £(£i,£2i*) and K,„(x) are defined as 

m,&,t) = - [(/ e u (*) - if (t)) c? + (i e 22 w - ifm el 1 



Vp, g (x) = <%[/ p , g (x) = //p// 2 cos(^(x - a))cos(/x 2 (y- c)), (p,g) £ Tmk- 



4.2.3. Approximations in 3D. In 3D case, again the function it(x, i n ) in 
also satisfies the Poisson equation in (|1.6[) which can be truncated on the computa- 
tional domain T> with homogeneous Dirichlet boundary conditions as 

(4.21) - V 2 u(x,t„) = |0(x,i„)| 2 , x£P; u(x, t n )|^ = 0. 

The above problem can be discretized by using a sine pseudospectral method in which 
the 0-modes are avoided. Denote the index set 

T M KL = {{p,q,r)\l<p<M-l, 1 < q < K - 1, 1 < r < L - 1} 

where M, K, L > are integers and define the functions 

Up, q , r (x) = sm(fi p (x - a))sin(/i 2 (y - c))sin(^(z - e)), (p,q,r) £ Tmkl, 
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where 

l4 = W(/ - e), l<r<L-l. 
Again, we take the (approximate) ansatz 

M-l K-l L-l 

(4.22) ti(x, t n ) = ^ ^ Kqr(tn) ^p,9,r(x), X = (x, 2) G X>, 

p— 1 q— 1 r— 1 

where u s pqr (t n ) is the sine transform of u(x, i„) corresponding to frequencies (p,q,r). 
Substituting (|4.22[) into the Poisson equation (|4.21[) and noticing the orthogonality of 
the sine functions, we obtain 



(4.23) K (t n ) = (p,q,r) GTmkl, 



where (\4> n \ 2 )pq r i- s the sine transform of |0(x,i„)| 2 corresponding to frequencies 

(p,q,r). 

Combining ()4.10j) . f|4. 22f) and ()4.23|) . we obtain an approximation of $(x, t) in the 
solution (|4.8[) via sine spectral method as 

M-lK-lL-l r 
p—1 q—1 r — 1 L 



(4-24) +e^l(X) + e^l(x) 



where the functions L(£i, £ 2 , 6, t), Vp,V,r(x), Vp ( , 2) , r (x) and V^, r (x) (for (p,q,r) G 
Tmkl) are defined as 

i(ei,6,6,t) = - [; e n W£ 2 +^ 22 we 2 2 +ewe 2 ] , 

^V.rC^) = d 5yU P ,qA*) = V v v\ cos(^(x - a)) cos(/x 2 (y - c)) sin(^(z- e)), 

^p (2 g ) ,r( x ) = fe^P.g.^x) = Mp^r cos(/i^(5 - a)) sin(/i 2 (y - c)) cos(^(z - e)), 

V p 3 q,r(*) = %?^p, 9 ,r(x) = /i 2 /i;! sin(/i£ (x - a)) cos^ 2 (y - c)) cos(/^ (z - e)). 

Remark 4.1. After obtaining the numerical solution 0(x,t) on the bounded 
computational domain T> , if it is needed to recover the original wave function ipfa, t) 
over a set of fixed grid points in the Cartesian coordinates x ; one can use the standard 
Fourier/sine interpolation operators from the discrete numerical solution cj)(5t,t) to 
construct an interpolation continuous function over T> \16\\41^ , which can be used to 
compute ^(x, t) over a set of fixed grid points in the Cartesian coordinates x for any 
fixed time t > 0. 

Remark 4.2. If the potential V(x) in {1.1$ is replaced by a time- dependent 
potential, e.g. V(x.,t), the rotating Lagrangian coordinates transformation and the 
numerical method are still valid provided that we replace W(x, t) in {3. 7| ) by W(x., t) = 
V(A(t)x, t) for x G M d andt> 0. 

5. Numerical results. In this section, we first test the accuracy of our numer- 
ical method, where throughout we apply the two-dimensional surface density model. 
Then study the dynamics of rotating dipolar BECs, including the center of mass, 
angular momentum expectation and condensate widths. In addition, the dynamics of 
vortex lattices in rotating dipolar BEC are presented. 
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Ax = 1/2 


Ax = 1/4 


Ax = 1/8 


Ax = 1/16 


/3 = 33.5914 


6.1569E-2 


1.7525E-4 


5.8652E-11 


<1E-11 


13 = 58.7849 


1.9746E-1 


2.3333E-3 


2.5738E-8 


2.6124E-11 


(3 = 92.3762 


4.8133E-1 


1.3385E-2 


1.6620E-6 


6.2264E-10 


j3 = 119.8488 


1.2984 


7.7206E-2 


9.5202E-5 


3.0974E-8 



Table 5.1: Spatial discretization errors \\<j){t) - (j) {Ax ' Ay ' At \t)\\ at time t = 1. 





At = 1/40 


At = 1/80 


At = 1/160 


At = 1/320 


At = 1/640 


/3 = 33.5914 


1.0434E-3 


2.6018E-4 


6.4992E-5 


1.6233E-5 


4.0456E-6 


/3 = 58.7849 


2.5241E-3 


6.2783E-4 


1.5674E-4 


3.9143E-5 


9.7550E-6 


13 = 92.3762 


4.9982E-3 


1.2380E-3 


3.0882E-4 


7.7108E-5 


1.9215E-5 


13 = 119.8488 


1.1417E-2 


2.7716E-3 


6.9009E-4 


1.7223E-4 


4.2915E-5 



Table 5.2: Temporal discretization errors \\4>(t) - ^ Aa: » A ^ At )(t)|| at time t = 1. 



5.1. Numerical accuracy. In order to test numerical accuracy, we consider a 
2D GPE (|4~l ) - (l4~2|) with the SDM long-range interaction ([LIB"]) and harmonic poten- 
tial ([1.19)1 . i.e. d = 2 in the GPE (|4.1[) . The other parameters are chosen as il = 0.4, 
lx = 7j/ = I, Tj = and dipole axis n = (0, 0, 1) T . The initial condition in (j4.3[) is 
taken as 

(5.1) 0o(x) = "174 e^ 2 ^, xeP, 

where wc perform our simulations on the bounded computational domain T> = [—16, 16] 2 
Denote (f> t - Ax ' A y' At > (t) as the numerical solution at time t obtained with the mesh size 
(Ax, Ay) and time step At. With a slight abuse of notation, wc let <f>(t) represent 
the numerical solution with very fine mesh size Ax = Ay = 1/64 and small time 
step At = 0.0001 and assume it to be a sufficiently good representation of the exact 
solution at time t. 

Tables I5.1H5.2I show the spatial and temporal errors of our numerical method for 
different j3 in the GPE (@TT]), where the errors are computed as ||0(t)-^ (A ^ Aff ' At) (i)||/2 
(with Ax = Ay) at time t = 1. To calculate the spatial errors in Table [5TT1 we always 
use a very small time step At = 0.0001 so that the errors from time discretization 
can be neglected compared to those from spatial discretization. Table 15.11 shows that 
the spatial accuracy of our method is of spectral order. In addition, the spatial errors 
increase with the nonlinearity coefficient (3 when the mesh size is kept constant. 

In Table 15.21 we always use mesh sizes Ai" = Ay = 1 /64 which are the same 
as those used in obtaining the 'exact' solution, so that one can regard the spatial 
discretization as 'exact' and the only errors are from time discretization. For different 
/3, Table [5~2l shows second order decrease of the temporal errors with respect to time- 
step size At. Similarly, for the same At, the temporal errors increase with f3. 

5.2. Dynamics of center of mass. In the following, we study the dynamics 
of the center of mass by directly simulating the GPE ([134} - (fl~T5)) in 2D with SDM 
long-range interaction (|1.13|) and harmonic potential ()1.19[) . To that end, we take 
d = 2, P = 30-\/l0/7r, rj = and dipole axis n = (1, 0, 0) T . The initial condition in 
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(|1.18j) is taken as 

(5.2) 0o (x) =aC(x-x°), with £(x) = (x + iy)e ~ <X ^ \ xeD, 

where the constant a is chosen to satisfy the normalization condition ||V'ol| 2 = !• 
Initially, we take x° = (1, 1) T . In our simulations, we use the computational domain 
T> = [—16, 16] 2 , the mesh size Ax = Ay = 1/16 and the time step size At = 0.0001. 

We consider the following two sets of trapping frequencies: (i) 7 X = "f y = 1, and 
(ii) "f x = 1, 7 y = 1.1. Figure [5TT1 shows the trajectory of the center of mass x c (£) in 
the original coordinates as well as the time evolution of its coordinates for different 
angular velocities f2, where "f x = "f y = 1. On the other hand, Figure IST21 presents 
the same quantities for y x = 1 and j y = 1.1. In addition, the numerical results are 
compared with analytical ones from solving the ODEs in (|2.14p - (|2.16|) . Figs. 15.1115.21 
show that if the external trap is symmetric, i.e. 7^ = 7^, the center of mass always 
moves within a bounded region which is symmetric with respect to the trap center 
(0,0) T . Furthermore, if the angular velocity f2 is rational, the movement is periodic 
with a period depending on both the angular velocity and the trapping frequencies. In 
contrast, when 7^ 7^ 7 y , the dynamics of the center of mass become more complicated. 
The simulation results in Figs. I5.1H5.2I arc consistent with those obtained by solving 
the ODE system in Lemma l2.3l for given Q, j x , and 7 y |51j and those numerical results 
reported in the literatures by other numerical methods [TDl[T2llT5] . 

On the other hand, we also study the dynamics of the center of mass x c (t) in the 
new coordinates. When 7^ = j y and D, arbitrary, the center of mass has a periodic 
motion on the straight line segment connecting — x° and x°. This is also true for 
x c (i) with f2 = (cf. Fig. 15.11) . However, the trajectories arc different for different f2 
if lx 7^ 7y • This observations agree with the results in Lemma 12.41 

In addition, our simulations show that the dynamics of the center of mass are 
independent of the interaction coefficients p and r], which is consistent with Lemma 
IO 

5.3. Dynamics of angular momentum expectation and condensate widths. 

To study the dynamics of the angular momentum expectation and condensate widths, 
we adapt the GPE ([iTT3)) - ([rT5]l in 2D with SDM long-range interaction ([TlB]) and 
harmonic potential (|1.19p . i.e. we take d — 2 and il = 0.7. Similarly, the initial 
condition in (| 1 . 18[) is taken as 

(5.3) ^o(x)=aC(x), xeP, 

where £(x) is defined in (|5.2p and a is a constant such that ||"0o|| 2 = 1- In our 
simulations, we consider the following four cases: 

(i) 7x = 7y = 1, p = 25^10^, V = 0, and n = (1, 0, 0) T ; 

(ii) j x = 7„ = 1, = 25^107^, r, = -15, and n = (1, 0, 0) T ; 

(iii) lx = ly = 1, p = 55^10/tt, t) = -15, and n = (0,0, 1) T ; 

(iv) lx = 1, j y = 1.1, p = 55v/l07^, V = -15, and n = (0, 0, 1) T . 

In Figure 15.31 we present the dynamics of the angular momentum expectation, 
energy and mass for each of the above four cases in the interval t £ [0, 15]. We see 
that if the external trap is radially symmetric in 2D, then the angular momentum 
expectation is conserved when either there is no dipolar interaction (Case (i)) or the 
dipolar axis is parallel to the z-axis (Case (iii)). Otherwise, the angular momentum 
expectation is not conserved. The above numerical observations are consistent with 



METHOD FOR ROTATING BOSE-EINSTEIN CONDENSATES 



19 




Fig. 5.1: Results for j x = j y = 1. Left: trajectory of the center of mass, x c (i) = 
(x c (t),y c (t)) T for < t < 100. Right: coordinates of the trajectory x c (i) (solid 
line: x c (t), dashed line: y c (t)) for different rotation speed fi, where the solid and 
dashed lines are obtained by directly simulating the GPE and '*' and 'o' represent 
the solutions to the ODEs in Lemma [2731 
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Fig. 5.2: Results for 7 X = l,7 y = 1.1. Left: trajectory of the center of mass, x c (t) = 
(x c (t),y c (t)) T for < t < 100. Right: coordinates of the trajectory x c (t) (solid 
line: x c (t), dashed line: y c (t)) for different rotation speed fi, where the solid and 
dashed lines are obtained by directly simulating the GPE and '*' and 'o' represent 
the solutions to the ODEs in Lemma [231 
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Case i 

- - Case ii 

■ — ■ Case iii 

■ ■ ■ ■ Case iv 



Energy: Case i 

— — Energy: Case ii 
■ — 1 - Energy: Case iii 

Energy: Case iv 

-# Mass: Case i-iv 



Fig. 5.3: Time evolution of the angular momentum expectation (left) and energy and 
mass (right) for Cases (i)-(iv) in section 5.3. 




the analytical results obtained in Lemma |2. II In addition, we find that our method 
conserves the energy and mass very well during the dynamics (cf. Fig. 15.31 right). 
Furthermore, from our additional numerical results not shown here for brevity, we 
observed that the angular momentum expectation is conserved in 3D for any initial 
data if the external trap is cylindrically symmetric and either there is no dipolar 
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interaction or the dipolar axis is parallel to the z-axis, which can also be justified 
mathematically. 

The dynamics of the condensate widths are presented in Figure 15.41 We find that 
S r (t) is periodic as long as the trapping frequencies satisfy ^ x = 7 y and the influence 
of the dipole axis vanishes, e.g. in the Case (i), which confirms the analytical results 
of Lemma 12.21 Furthermore, from our additional numerical results not shown here 
for brevity, we observed that S r (t) is periodic and 8 x (t) = 5 y (t) = \8 r {t) if -q = for 
any initial data or n = (0,0, 1) T for radially symmetric or central vortex-type initial 
data. 



t = t = 1 .5 t = 3 




Fig. 5.5: Contour plots of the density function ^(x, t)\ 2 for dynamics of a vortex 
lattice in a rotating BEC (Case (i)). Domain displayed: (x,y) € [—13, 13] 2 . 



5.4. Dynamics of quantized vortex lattices. In the following, we apply our 
numerical method to study the dynamics of quantized vortex lattices in rotating 
dipolar BECs. Again, we adapt the CPE ([TTTl)) - (jTTT5]) in 2D with SDM long-range 
interaction (|1.13[) and harmonic potential (|1.19l) . i.e. we choose d = 2, /3 = 1000 and 
H = 0.9. The initial datum in (|1.18[) is chosen as a stationary vortex lattice which is 
computed numerically by using the method in |49[I50| with the above parameters and 
"fx — ly = 1, V = 0, i.e. no long-range dipolc-dipole interaction initially. Then the 
dynamics of vortex lattices are studied in two cases: 

(i) perturb the external potential by setting ^ x = 1.05 and j y = 0.95 at t = 0; 

(ii) turn on the dipolar interactions by setting rj = —600 and dipolar axis n = 
(1,0, 0) T at time t = 0. 

In our simulations, we use T> = [—16, 16] 2 , Ax = Ay = 1/16 and At = 0.0001. Figures 
15.51 15. 6l show the contour plots of the density function ^(xji)! 2 at different time steps 
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for Cases (i) and (ii), respectively, where the wave function tp(x,t) is obtained from 
</>(x, t) by using interpolation via sine basis (see Remark 14.11) . We see that during 
the dynamics, the number of vortices is conserved in both cases. The lattices rotate 
to form different patterns because of the anisotropic external potential and dipolar 
interaction in Cases (i) and (ii), respectively In addition, the results in Case (i) are 
similar to those obtained in |10j . where a spectral type method in polar coordinates 
was used to simulate the dynamics of vortex lattices. 



t = t = 1 .5 t = 3.5 




Fig. 5.6: Contour plots of the density function |?/>(x, t)\ 2 for dynamics of a vortex 
lattice in a rotating dipolar BEC (Case (ii)). Domain displayed: (x, y) € [—10, 10] 2 . 



6. Conclusions. Wc proposed a simple and efficient numerical method to simu- 
late the dynamics of rotating dipolar Bosc-Einstcin condensation (BEC) whose prop- 
erties are described by the Gross-Pitaevskii equation (CPE) with both the angu- 
lar rotation term and the long-range dipolc-dipolc interaction. First, by decoupling 
the short-range and long-range interactions, we reformulate the GPE as a Gross- 
Pitaevskii- (fractional) Poisson system. Then wc eliminate the angular rotation term 
from the GPE using a rotating Lagrangian coordinate transformation, which makes 
it possible to design a simple and efficient numerical method. In the new rotat- 
ing Lagrangian coordinates, we presented a numerical method which combines the 
time-splitting techniques with Fourier/sine pseudospcctral approximation to simulate 
the dynamics of rotating dipolar BECs. The numerical methods is explicit, uncon- 
ditional stable, spectral accurate in space and second order accurate in time, and 
conserves the mass in the discretized level. The memory cost is 0(MK) in 2D and 
O(MKL) in 3D, and the computational cost per time step is O (MKln(MK)) in 
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2D and O (MKL\n(MKL)) in 3D. More specifically, the method is very easy to be 
implemented via FFT or DST. We then numerically examine the conservation of the 
angular momentum expectation and study the dynamics of condensate widths and 
center of mass for different angular velocities. In addition, the dynamics of vortex 
lattice in rotating dipolar BEC are investigated. Numerical studies show that our 
method is very effective in simulating the dynamics of rotating dipolar BECs. 
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